*==========================================================================================================================================================
*
*	SINGLE EQUATION REGRESSIONS
*
*==========================================================================================================================================================

clear all

set more off

qui{

*-------------------------------------------------------------------------------------------------------------
*	SETTINGS
*-------------------------------------------------------------------------------------------------------------

* Variables
local vars "gdp_r_adj_gr cpi_sil_gr log_R_esp_er_hp_gr itsN wage_sil_gr " 
local lava ""Real output" "Consumer prices" "Lending rate" "Money inflow (% of stock)" "Nominal wage"" 
local yti ""Percent" "Percent" "Percentage pts" "Percentage pts" "Percent"" 


* Lags
local var_lags "4 4 4 4 4" /*number of lags of dependend variable*/
local shock_lags "4 4 4 4 4" /*number of shock lags*/

* Bootstrap runs
local B 500 //500

* irf horizon
local H 5
local F 5

* maritime disaster IV
local dis "perclossN"

*-------------------------------------------------------------------------------------------------------------

local vecsize : word count `vars'
forvalues s = 1/`vecsize' { /*loop over indicators*/
local v : word `s' of `vars'

local progress0 = "" /*display progress report*/
n display "`progress0'"
local progress1 = "`v'"
n display "`progress1'"


* Number of lags
local L : word `s' of `var_lags' /* number of lagged dependend variable */
local L_shock : word `s' of `shock_lags' /* number of lagged shock variable */


***_________ CONTROLS
local indep	" gdp_r_adj_gr cpi_sil_gr wage_sil_gr itsN_gr"

forvalues i=0/`L'{
	local indep_`i' = subinstr("`indep'"," "," l`i'_",.)
	}
	
local rhs ""
forvalues i = 0/`L'{
	local rhs "`rhs'`indep_`i''"
	}
	
	local not l0_`v'_gr l0_`v' l0_itsN_gr 
	local rhs: list rhs- not
	*/


*_________ GET DATA ____________
use "data/working files/prepared.dta", clear /*read in data*/


*____ DECLARE TIME SERIES _____
tsset year, yearly



*_______________________pVAR-Analysis__________________________

*specification
local spec l(0/`L_shock').perclossN l(-`F'/-1).perclossN l(-`F'/`L_shock').perclossN_vellon l(-`F'/`L_shock').percloss_salvN l(-`F'/`L_shock').percloss_pirN l(0/`L_shock').temperature l(0/`L_shock').warbeu l(0/`L_shock').warbciv 

if(`L'>0){
reg `v' l(1/`L').`v' `spec'
	}
else{ /*no lagged dependent regressors*/
reg `v' `spec'
	}
	
		gen d = 1 if e(sample)==1
		sum d
		mat N = r(N)
			
mat est_1 = e(b)
predict res_`v'_1, residuals
	
* Bring into canonical form: zeta_t = F*zeta_(t-1) + G*S_(t)

if(`L'>=1){ /*only if lagged dependent variables are contained*/

*matrix F
	if(`L'>1){
		matrix h1 = J(`L'-1,1,1)
		matrix F = diag(h1)
		matrix h2 = J(`L'-1,1,0)
		matrix F = F,h2
		matrix F =est_1[1, 1..`L']\F
		}
	if(`L'==1){ /*1x1 matrix*/
		matrix F = est_1[1,1]
		}
	}

*matrix G
local LL = `L'+1
local LLL = `L' + `L_shock' +1

if(`L'>1){
	matrix G = est_1[1, `LL'..`LLL']
	matrix h = J(`L'-1,`L_shock'+1,0)
	matrix G = G\h
	}
else{ /*i.e. one or no lagged dependent variable*/
	matrix G = est_1[1, `LL'..`LLL']
	}

*_____________________________IMPULSE RESPONSE FUNCTIONS______________________

if(`L'>=1){ /*only if lagged dependent variables are contained*/

* F-exponential matrices
local expF = "F"
matrix F0 = 1
matrix F1 = F

forvalues x=2/`H'{
	local expF_`x' = subinstr(" F", " ", "*", .)
	local expF = subinstr("`expF' ", " ", "`expF_`x''", 1)
	matrix F`x' = `expF'
	}
	
forvalues x=0/`H'{ /*obtain "FG" matrices*/
	matrix F`x'G = F`x'*G
	}

* marginal effect-matrix Psi: Psi*S = IRF
matrix Psi = J(`H'+1,`H'+1,0)

local h = `H'+1
forvalues i = 1/`h'{ /*rows*/
	forvalues j = 1/`h'{ /*columns*/
		local ij = `i'-`j'
		if(`j'<=`i' & `j'<=`L_shock'+1){ /*below and on the diagonal - entries*/ /*j<than number of lags of exogenous shock variable*/
			mat Psi[`i',`j'] = F`ij'G[1,`j']
			}
		}
	}
}
	
else{ /*i.e. no lagged dependent variable*/

matrix Psi = J(`H'+1,`H'+1,0)

local h = `H'+1
forvalues i = 1/`h'{ /*rows*/
	forvalues j = 1/1{ /*columns*/
		if(G[1,`i']!=.){
			mat Psi[`i',`j'] = G[1,`i']
			}
		else{
			mat Psi[`i',`j'] = 0
			}
		}
	}	
}

* shock matrices: S
forvalues x=1/`h'{
	matrix S`x' = J(`h',1,0)
	local y = `x'-1
	if(`x'>1){
		matrix S`x'[`y',1] = 0 /*set previous entry to 0*/
		}
	matrix S`x'[`x',1] = 1 /*set current entry to 1*/
	}

* IRFs
matrix irf = J(`H'+1,1,0)
	
forvalues x=1/6{
	//local x2 = `x'+1
		
	matrix irf_`x' = Psi'*S`x'
	
	mata: irf_`x' = st_matrix("irf_`x'")
	mata: irf_`x' = sum(irf_`x')
	mata: st_matrix("irf_`x'", irf_`x')
	
	matrix irf[`x',1] = irf_`x'
	}

* cumulate and save as variable
svmat irf, names(irf_`v')
if (`"`v'"'=="itsN") {
	gen irf_`v' = irf_`v'1 if _n<=`H'+1
	}
else {
	gen irf_`v' = sum(irf_`v'1) if _n<=`H'+1
	}

svmat N, names(N_`v')
ren N_`v'1 N_`v'

drop irf_`v'1










*___________________BOOTSTRAP - CONFIDENCE INTERVALS__________________

* resampling
mkmat res_`v'_1, matrix(res_`v'_1) nomissing
svmat res_`v'_1, name(res_`v') /*residual variable without missings*/

local a = 1
local b = rowsof(res_`v'_1)

** draw residuals with replacement
forvalues x=1/`B' {
	set seed `x'
	generate rnd_`x' = floor((`b'-`a'+1)*runiform()+`a') /*draw random integer for residual selection*/
	generate error_`x' = res_`v'1[rnd_`x'[_n]]
	}

** fitted values
if(`L'>0){
reg `v' l(1/`L').`v' `spec'
	}
else{ /*no lagged dependent regressors*/
reg `v' `spec'
	}
				
predict fit

** new samples: fitted + resampled error
forvalues x=1/`B' {
	gen `v'_`x' = fit + error_`x'
	}

* IRFs
local count = " 0"

forvalues b = 1/`B' {

local count = `count'+1
local count2 = " `count'"
local progress2 = subinstr("`count2'"," ","Bootstrap ", .)
n display "`progress2'"
	
if(`L'>0){
reg `v'_`b' l(1/`L').`v' `spec'
	}
else{ /*no lagged dependent regressors*/
reg `v'_`b' `spec'
	}
		
mat est_1 = e(b)

* Bring into canonical form: zeta_t = F*zeta_(t-1) + G^np*S_(t-1) + G^p*S_(t-1)

if(`L'>=1){ /*only if lagged dependent variables are contained*/

*matrix F
	if(`L'>1){
		matrix h1 = J(`L'-1,1,1)
		matrix F = diag(h1)
		matrix h2 = J(`L'-1,1,0)
		matrix F = F,h2
		matrix F =est_1[1, 1..`L']\F
		}
	if(`L'==1){ /*1x1 matrix*/
		matrix F = est_1[1,1]
		}
	}

*matrix G
local LL = `L'+1
local LLL = `L' + `L_shock' +1

if(`L'>1){
	matrix G = est_1[1, `LL'..`LLL']
	matrix h = J(`L'-1,`L_shock'+1,0)
	matrix G = G\h
	}
else{ /*i.e. one or no lagged dependent variable*/
	matrix G = est_1[1, `LL'..`LLL']
	}
		
*_____________________________IMPULSE RESPONSE FUNCTIONS______________________

if(`L'>=1){ /*only if lagged dependent variables are contained*/

* F-exponential matrices
local expF = "F"
matrix F0 = 1
matrix F1 = F

forvalues x=2/`H'{
	local expF_`x' = subinstr(" F", " ", "*", .)
	local expF = subinstr("`expF' ", " ", "`expF_`x''", 1)
	matrix F`x' = `expF'
	}
	
forvalues x=0/`H'{ /*obtain "FG" matrices*/
	matrix F`x'G = F`x'*G
	}

* marginal effect-matrix Psi: Psi*S = IRF
matrix Psi = J(`H'+1,`H'+1,0)

local h = `H'+1
forvalues i = 1/`h'{ /*rows*/
	forvalues j = 1/`h'{ /*columns*/
		local ij = `i'-`j'
		if(`j'<=`i' & `j'<=`L_shock'+1){ /*below and on the diagonal - entries*/ /*j<than number of lags of exogenous shock variable*/
			mat Psi[`i',`j'] = F`ij'G[1,`j']
			}
		}
	}
	}
	
else{ /*i.e. no lagged dependent variable*/

matrix Psi = J(`H'+1,`H'+1,0)

local h = `H'+1
forvalues i = 1/`h'{ /*rows*/
	forvalues j = 1/1{ /*columns*/
		if(G[1,`i']!=.){
			mat Psi[`i',`j'] = G[1,`i']
			}
		else{
			mat Psi[`i',`j'] = 0
			}
		}
	}	
	}

* shock matrices: S
forvalues x=1/`h'{
	matrix S`x' = J(`h',1,0)
	local y = `x'-1
	if(`x'>1){
		matrix S`x'[`y',1] = 0 /*set previous entry to 0*/
		}
	matrix S`x'[`x',1] = 1 /*set current entry to 1*/
	}

* IRFs
matrix irf = J(`H'+1,1,0)
	
local HH = `H'+1
forvalues x=1/`HH'{
	//local x2 = `x'+1
		
	matrix irf_`x' = Psi'*S`x'
	
	mata: irf_`x' = st_matrix("irf_`x'")
	mata: irf_`x' = sum(irf_`x')
	mata: st_matrix("irf_`x'", irf_`x')
	
	matrix irf[`x',1] = irf_`x'
	}

	
* cumulate and save as variable
svmat irf, names(irf_`v'_`b'_)
if (`"`v'"'=="itsN") {
	gen irf_`v'_cum`b' = irf_`v'_`b'_1 if _n<=`H'+1
	}
else {
	gen irf_`v'_cum`b' = sum(irf_`v'_`b'_1) if _n<=`H'+1
	}

drop `v'_`b' error_`b'

}

** 95 Percentiles of IRFs
***lower bound
egen ci_lo_`v' = rowpctile(irf_`v'_cum*), p(15.9)
egen ci_lo90_`v' = rowpctile(irf_`v'_cum*), p(5)

***upper bound
egen ci_up_`v' = rowpctile(irf_`v'_cum*), p(84.1)
egen ci_up90_`v' = rowpctile(irf_`v'_cum*), p(95)
	
drop irf_`v'_cum* irf_`v'_*_1
	
* plot prerequisites
gen time = _n-1 if(_n<=`H'+1) 

egen N = min(N_`v')
drop N_`v'
ren N N_`v'	
drop if time == . 
gen zero = 0 if time!=.

local title : word `s' of `lava'
local ytitle : word `s' of `yti'
local lsize "medlarge"
local tsize "medlarge"

twoway  (rarea ci_lo90_`v' ci_up90_`v' time, fcolor(gs15) lcolor(gs15)) ///
		(rarea ci_lo_`v' ci_up_`v' time, fcolor(gs12) lcolor(gs12)) ///
		(line irf_`v' time, lcolor(black) lpattern(solid) lwidth(medthick)) ///
		(line zero time, lcolor(black) lwidth(medium)), title("") xtitle("Year", size(`lsize')) ytitle(`ytitle', height(2) size(`lsize')) xlabel(,labsize(`lsize')) ylabel(, format(%9.1f) nogrid angle(0) labsize(`lsize')) ///
		legend(off) title(`title', size(`tsize')) ///
		graphregion(color(white)) plotregion(color(white)) name(g`s', replace) nodraw
			
		

	} /*variables*/
	

graph combine g4 g1 g3 g2 g5, name(SE_presi, replace) cols(3) imargin(medsmall) plotregion(color(white)) graphregion(color(white)) 
graph display SE_presi, ysize(4) xsize(6)
graph export "results/FigureB9.pdf", replace 

	
}
